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We  report  a  method  that  enables  long-time  molecular  dynamics  (MD)  simulations  of  shock  wave  loading. 

The  goal  is  to  mitigate  the  severe  interference  effects  that  arise  at  interfaces  or  free  boundaries  when  using 
standard  nonequilibrium  MD  shock  wave  approaches.  The  essence  of  the  method  is  to  capture  between  two 
fixed  pistons  the  material  state  at  the  precise  instant  in  time  when  the  shock  front,  initiated  by  a  piston  with 
velocity  Up  at  one  end  of  the  target  sample,  traverses  the  contiguous  boundary  between  the  target  and  a  second, 
stationary  piston  located  at  the  opposite  end  of  the  sample,  at  which  point  the  second  piston  is  also  assigned 
velocity  Up  and  the  simulation  is  continued.  Thus,  the  target  material  is  captured  in  the  energy-volume  Hugo- 
niot  state  resulting  from  the  initial  shock  wave,  and  can  be  propagated  forward  in  time  to  monitor  any 
subsequent  chemistry,  plastic  deformation,  or  other  time-dependent  phenomena  compatible  with  the  spatial 
scale  of  the  simulation.  For  demonstration  purposes,  we  apply  the  method  to  shock-induced  chemistry  in 
methane  based  on  the  adaptive  intermolecular  reactive  empirical  bond  order  force  field  [S.  J.  Stuart  et  al,  J. 

Chem.  Phys.  112,  6472  (2000)]. 

DOI:  10.1103/PhysRevB.76.224108  PACS  number(s):  62.50.+p,  82.20.Wt,  82.40.Fp 


I.  INTRODUCTION 

Dynamic  material  response  to  shock  wave  loading  has 
been  studied  for  decades  for  both  practical  and  fundamental 
reasons.^  Molecular  dynamics  (MD)  is  the  most  widely  used 
method  for  theoretical  studies  of  physical  and  chemical  pro¬ 
cesses  in  condensed  materials  on  submicron  scales.  Since  the 
1970s,  nonequilibrium  molecular  dynamics  (NEMD)  has 
been  applied  to  studies  of  shock-induced  phenomena  such  as 
defect  generation,^  phase  transitions,^  and  chemistry.'^’^  Al¬ 
though  current  computer  capabilities  allow  shock  wave 
simulations  for  systems  containing  even  billions  of  atoms 
interacting  via  comparatively  simple  potentials  (for  instance, 
the  embedded- atom  model), ^  it  is  computationally  expensive 
to  use  NEMD  to  study  shock  waves  in  complicated,  poly¬ 
atomic  molecular  crystals  characterized  by  many-body  inter¬ 
actions  and  electrostatic  contributions  to  the  potential 
energy.^  This  is  particularly  true  since  many  of  the  phenom¬ 
ena  of  interest  occur  on  time  scales  significantly  larger  than 
those  required  for  the  passage  of  a  shock  wave  through  a 
typical  MD  simulation  cell. 

Various  equilibrium  MD  methods,  sometimes  called 
Hugoniostat  methods,  have  been  developed  in  the  past  de¬ 
cade  to  reproduce  the  final  states  of  shocked  materials  by 
way  of  extended  equations  of  motion  that  act  on  the  system 
to  drive  it  toward  a  prescribed  Hugoniot  state. (The 
Hugoniot  is  the  locus  of  states  accessible  by  shock  wave 
loading.)  This  allows  long-time  sampling  of  shock  states 
without  the  need  to  simulate  a  system  that  is  large  enough  to 
sustain  the  compressed  state  for  a  long  enough  time  to  arrive 
at  those  states.  In  cases  where  the  shock  wave  dynamics  is  of 
interest,  however,  the  Hugoniostat  methods  are  of  limited 
use.  Zhakhovskii  et  al}^  developed  a  “moving  window” 
method  that  allows  the  study  of  the  detailed  dynamics  in  the 


vicinity  of  a  shock  front  by  systematically  adding  new  ma¬ 
terial  to  the  unshocked  end  of  the  simulation  cell  and  remov¬ 
ing  shocked  material  from  the  other  end  so  that  the  simula¬ 
tion  size  remains  constant  and  the  simulation  frame  remains 
centered  on  the  shock  front.  However,  even  this  method  be¬ 
comes  difficult  to  apply  to  simulations  with  nonprompt 
chemistry  or  mechanical  deformation,  since  adding  and  re¬ 
moving  material  from  the  simulation  will  certainly  introduce 
artifacts  for  practically  accessible  simulation  domains.  Most 
recently,  Zhao  et  al  developed  a  NEMD  method  to  study 
shock-induced  alloying  reactions  in  Ni/Al  nanolaminates.  In 
their  method,  shock  waves  are  first  generated  by  colliding 
two  identical  slabs  with  equal  but  opposite  center-of-mass 
velocities.  The  periodic  boundaries  along  the  shock  direction 
shrink  consistently  with  the  mass  velocity.  Once  the  shock 
waves  reach  the  periodic  boundaries,  however,  the  cell  pa¬ 
rameters  are  fixed,  and  the  simulation  is  continued  in  the 
usual  three-dimensional  periodic,  constant  energy-constant 
volume  (NVE)  ensemble  to  allow  the  alloying  reactions  to 
approach  equilibrium. 

Here,  we  report  a  NEMD  method  that  allows  the  dynam¬ 
ics  behind  a  shock  wave  to  evolve  with  minimal  interference 
from  the  free  surface  of  the  simulation  cell.  We  demonstrate 
the  method  in  a  study  of  shock-induced  chemistry  in 
condensed-phase  methane  based  on  the  adaptive  intermo¬ 
lecular  reactive  empirical  bond  order  (AIREBO)  force 
field. 


II.  COMPUTATIONAL  DETAILS 

A.  Absorbing  boundary  condition 

A  schematic  diagram  of  how  the  absorbing  boundary  con¬ 
dition  is  applied  in  a  simulation  of  a  shock  wave  is  presented 
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FIG.  1.  (Color  online)  Schematic  illustration  of  how  the  shock 
absorbing  boundary  condition  is  applied,  (a)  A  shock  wave  is  gen¬ 
erated  in  the  sample  by  driving  a  rigid  piston  into  it  with  constant 
velocity  Up\  a  second  piston  is  contiguous  to  and  equilibrated  with 
the  material  on  the  opposite  end  of  the  simulation  cell,  (b)  The 
shock  wave  reaches  the  second  piston,  (c)  The  second  piston  begins 
to  move  with  the  same  fixed  velocity  as  the  first. 

in  Fig.  1.^^  By  moving  a  rigid  piston  (piston- 1)  at  constant 
velocity  Up,  a  shock  wave  is  generated  that  moves  with  ve¬ 
locity  Us  through  the  target  material  [see  Fig.  1(a)].  On  the 
opposite  end  of  the  simulation  cell  is  a  second  rigid  piston 
(piston-2),  contiguous  with  the  target  and  assigned  zero  ve¬ 
locity.  When  the  shock  wave  reaches  piston-2  [Fig.  1(b)],  it 
is  instantaneously  assigned  the  same  constant  velocity  as 
piston- 1  {up)  [see  Fig.  1(c)].  From  that  point  onward,  the 
simulation  is  microcanonical  and  maintains  the  initial  Hugo- 
niot  state  associated  with  the  passage  of  the  shock  front. 
Chemical  reactions  or  other  dynamical  processes  can  be  fol¬ 
lowed  until  they  reach  equilibrium.  These  processes  can  re¬ 
sult  in  significant  changes  in  temperature,  pressure,  and  com¬ 
position  in  the  confined  region. 

A  critical  issue  is  exactly  when  and  how  to  apply  the 
velocity  Up  to  the  second  piston.  Although  various  criteria 
can  be  imagined,  the  initial  transfer  of  internal  energy  is  a 
reasonable  one  for  defining  the  instant  of  shock  wave  pas¬ 
sage  across  a  given  dividing  surface  in  configuration  space 
since  shock  wave  propagation  is  essentially  energy  transfer 
from  the  moving  piston  to  lattice  degrees  of  freedom  in  the 
target  material.  We  show  in  Fig.  2  how  we  use  the  internal 
energy  to  determine  the  time  at  which  the  second  piston  be¬ 
gins  to  move.  In  this  case,  the  piston  velocity  Up  =  3.0  km/s 
is  below  the  threshold  for  shock-induced  chemistry  in  meth¬ 
ane.  The  simulation  cell  is  arbitrarily  divided  into  80  bins, 
each  one  unit  cell  wide;  these  bins  deform  affinely  as  the 
shock  wave  compresses  the  sample.  Figure  2(a)  contains  the 
internal  energy  profile  in  the  material  along  the  shock  direc¬ 
tion  at  t=5A5  ps,  which  is  just  before  the  internal  energy  of 
the  bin  immediately  adjacent  to  piston-2  begins  to  rise  rap¬ 
idly  [denoted  as  bin  80,  corresponding  to  zero  displacement 
along  the  abscissa  in  Fig.  2(a);  the  internal  energy  of  this  bin 
is  denoted  as  ego  in  Fig.  2(b)].  The  internal  energy  profile  of 
the  shocked  material  in  Fig.  2(a)  is  fitted  to  a  straight  line, 
which  is  extrapolated  to  predict  the  internal  energy  SgQ  at  the 
moment  when  the  shock  front  just  reaches  its  outer  edge 
(that  is,  the  contiguous  boundary  with  the  second  piston). 
Monitoring  the  internal  energy  SgQ  as  the  trajectory  contin¬ 
ues,  we  find  that  the  extrapolated  value  is  achieved  at  t 
=  5.57  ps  [see  Fig.  2(b)];  the  dashed  horizontal  line  in  Fig. 
2(b)  shows  the  extrapolated  value  for  the  internal  energy 
corresponding  to  the  moment  when  the  shock  wave  reaches 
piston-2.  It  is  at  this  time  that  the  second  piston  begins  to 
move  with  the  same  fixed  velocity  Up  as  the  first  one. 


FIG.  2.  Summary  of  details  for  determining  when  the  second 
piston  begins  to  move,  (a)  Internal  energy  profile  along  the  shock 
direction;  (b)  time  evolution  of  internal  energy  in  the  material  sub¬ 
volume  (bin)  closest  to  the  second  piston  [zero  displacement  along 
the  abscissa  in  (a)  corresponds  to  internal  energy  sgQ  in  (b)].  The 
time  snapshot  in  (a)  is  for  t=5.45  ps;  the  line  is  the  linear  fit  for  the 
shocked  material.  The  dashed  horizontal  line  in  (b)  shows  the  ex¬ 
trapolated  value  for  the  internal  energy,  corresponding  to  the  mo¬ 
ment  when  the  shock  wave  reaches  piston-2.  When  the  internal 
energy  sgQ  reaches  this  value  {t=5.51  ps),  the  second  piston  is  as¬ 
signed  a  constant  velocity  Up. 

There  is  a  significant  disparity  in  the  time  scales  for  shock 
wave  traversal  of  the  sample  in  most  NEMD  simulations  and 
subsequent  establishment  of  chemical  or  thermo-mechanical 
equilibrium.  In  the  simplest  sense,  the  maximum  time  acces¬ 
sible  to  the  former  is  the  shock  transit  time  across  the 
sample,  t^ax= I  sample  ^  However,  this  only  applies  to  mate¬ 
rial  in  the  immediate  vicinity  of  the  first  piston;  material  on 
the  free  boundary  is  under  compression  for  essentially  zero 
time.  (The  time  required  for  a  backscattered  wave  to  re¬ 
traverse  the  compressed  system  sets  the  true  upper  limit  on 
the  time  that  any  region  in  such  a  simulation  can  be  sustained 
in  the  shocked  state.)  The  proposed  absorbing  boundary  con¬ 
dition  described  here  minimizes  in  a  practical  way  the  effects 
of  wave  reflection  from  a  free  surface,  effectively  providing 
a  near-perfect  impedance  match^^  between  the  target  material 
and  the  two  pistons,  and  thus,  allows  simulation  of  the 
sample  in  a  shock-compressed  state  for  an  interval  of  time 
whose  limit  is  determined  by  the  stability  of  the  numerical 
integration  scheme. 

B.  Model  system  and  details  of  simulation  procedure 

For  demonstration  purposes,  we  have  chosen  to  study 
shock-induced  chemistry  in  methane  as  predicted  by  the 
AIREBO  potential  due  to  Stuart  et  aO^  AIREBO  is  an  ex¬ 
tension  of  the  reactive  empirical  bond  order  potential. 
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AIREBO  was  designed  to  optimally  describe  liquid- state  hy¬ 
drocarbon  properties  at  ambient  pressure,  while  perturbing  as 
little  as  possible  predictions  for  the  solid-state  polymorphs  of 
pure  carbon.  AIREBO  has  been  used  previously  in  MD  stud¬ 
ies  of  thermal  dissociation  in  methane,  ethylene,  and  ben¬ 
zene,  and  of  shock-induced  chemistry  in  solid  ace- 
tylene,^^’^^  ethylene,^^  methane,^^’^^  and  anthracene.^^  The 
MD  simulations  presented  here  were  performed  using  a  com¬ 
puter  code  developed  for  AIREBO  by  Stuart  et  al  Trajecto¬ 
ries  were  integrated  using  the  velocity  Verlet  algorithm,  with 
step  sizes  in  the  interval  from  0.1  to  0.25  fs  depending  on 
the  temperature  and  pressure. 

The  simulation  cell  for  the  shock  wave  simulations  con¬ 
sists  of  82  X  3  X  3  unit  cells  of  methane  phase  I,  a  face- 
centered-cubic  (fee)  unit  cell  with  rotating  molecules  at  the 
lattice  sites. Initially,  all  atoms  were  positioned  on  the  per¬ 
fect  fee  lattice  with  no  orientational  disorder.  Next,  the  sys¬ 
tem  was  equilibrated  in  the  NVT  ensemble,  with  periodic 
boundary  conditions  applied  in  the  directions  transverse  to 
the  direction  of  subsequent  shock  loading,  and  cell  param¬ 
eters  adjusted  to  yield  zero  pressure  at  a  given  temperature;  a 
gap  of  11.722  A  (two  unit  cell  widths)  was  introduced  be¬ 
tween  the  first  and  second  layers  of  unit  cells  along  the  shock 
direction.  Atoms  in  the  first  and  last  layers  of  unit  cells  along 
the  shock  direction,  which  comprise  piston- 1  and  piston-2, 
respectively,  were  held  fixed  during  this  equilibration  period, 
which  was  continued  until  the  system  reached  steady- 
fluctuating  values  about  the  prescribed  temperature.  Shock 
waves  were  generated  by  driving  piston- 1  into  the  target 
sample  with  constant  velocity  Up.  Piston-2,  on  the  opposite 
end  of  the  sample,  was  initially  assigned  zero  velocity.  Both 
pistons  were  perfectly  rigid. 

III.  RESULTS  AND  DISCUSSION 

A.  Comparison  of  adaptive  intermolecular  reactive  empirical 
bond  order  predictions  to  experimental  results 

We  provide  details  of  the  validation  of  the  absorbing 
boundary  condition  in  Sec.  Ill  B.  In  this  section,  we  present 
a  comparison  between  AIREBO  predictions  obtained  using 
the  absorbing  boundary  condition  and  experimental  shock 
wave  data^^’^^  for  liquid  methane  at  initial  temperature 
111  K  and  density  0.42  g/cm^;  the  results  are  summarized  in 
Eig.  3.  The  shock  temperature  predicted  for  AIREBO  is  in 
good  agreement  with  experiment  [Eig.  3(a)].  Eigure  3(b)  in¬ 
dicates  that  the  AIREBO  potential  overestimates  the  shock 
speed  by  —25%  for  Up=^3  km/s,  whereas  the  compression 
ratio  p/po  in  Eig.  3(c)  is  underestimated  by  —30%  at  that 
same  piston  velocity.  Einally,  the  component  of  stress  along 
the  shock  direction  (hereafter  referred  to  as  shock  pres¬ 
sure)  behind  the  shock  front  is  overestimated  by  about  25% 
at  the  highest  piston  velocity  considered  [see  Eig.  3(d)];  this 
discrepancy  tends  to  zero  with  decreasing  piston  velocity. 
Since  temperature  is  the  dominant  thermodynamic  variable 
for  chemistry,  the  agreement  in  Eig.  3(a)  suggests  that  shock- 
induced  reactions  predicted  for  methane  by  AIREBO  may  be 
reasonable,  and  are  certainly  sufficient  for  the  present  goal  of 
methods  development. 


FIG.  3.  Shock  strength  dependence  of  (a)  temperature,  (b)  shock 
speed,  (c)  compression  ratio,  and  (d)  shock  pressure  in  liquid  meth¬ 
ane.  The  initial  temperature  and  density  are  111  K  and  0.42  g/em^, 
respectively.  Solid  symbols,  simulation  predictions,  open  symbols, 
experiment.  Data  for  pressure,  eompression  ratio,  and  shock  veloc¬ 
ity  are  from  Nellis  et  al  (Ref.  23);  temperature  measurement  is 
from  Radousky  et  al  (Ref  24). 

B.  Validation  of  the  absorbing  boundary  condition 

Figure  4  contains  plots  of  density,  material  velocity,  inter¬ 
nal  energy,  and  temperature  profiles  for  subvolumes  of  the 
simulation  cell  along  the  shock  direction,  at  times  before  and 
after  the  second  piston  begins  to  move.  The  results,  eight 
snapshots  in  time,  are  for  a  i/^  =  3.0  km/s  shock  in  crystalline 
methane  equilibrated  at  50  K  and  zero  pressure.  Traces  for 
successive  times  are  shifted  vertically  along  the  ordinate  as 
an  aid  to  the  eye.  In  this  case,  the  shock  wave  reaches  the 
second  piston  at  t=  —5.57  ps  (bold  trace  in  Fig.  4),  at  which 
time  piston-2  is  assigned  a  velocity  Up.  Whereas  this  would 
be  the  maximum  time  (or,  at  best,  half  the  maximum  time) 
accessible  by  other  NEMD  shock  methods,  with  the  excep¬ 
tion  of  that  of  Zhao  et  al,^^  in  the  present  case,  the  simula¬ 
tion  was  continued  for  an  additional  30  ps.  One  can  see  from 
Fig.  4  that  spatial  distributions  of  the  density,  local  mass 
velocity,  and  temperature  are  essentially  constant  over  the 
entire  time  interval  after  the  second  piston  starts  to  move. 
The  profiles  of  internal  energy — which  are  the  basis  for  de¬ 
termining  when  the  shock  wave  traverses  the  sample  bound¬ 
ary  into  the  second  piston — have  a  negative  slope  across  the 
simulation  cell  during  shock  wave  passage.  In  this  case,  the 
slope  is  approximately  preserved  immediately  after  the  sec¬ 
ond  piston  begins  to  move,  but  disappears  within  5  ps.  There 
is  no  evidence  in  any  of  the  results  shown  in  Fig.  4  for 
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Distance  from  piston-2  (A)  Distance  from  piston-2  (A) 


EIG.  4.  Spatial  profiles  of  (a)  density,  (b)  local  velocity,  (c) 
internal  energy,  and  (d)  temperature  along  the  shock  direction,  for 
eight  time  snapshots  before  and  after  the  absorbing  boundary  con¬ 
dition  is  applied  to  shocked  solid  methane.  The  piston  velocity  is 
3  km/  s  and  the  initial  temperature  is  50  K.  The  absorbing  boundary 
condition  is  applied  at  t=5.51  ps  (bold  trace  in  the  figures).  Traces 
for  successive  snapshots  are  shifted  vertically  for  clarity. 


significant  reflections  or  buildup  of  energy  at  either  piston- 
sample  boundary. 

In  Fig.  5,  we  compare  results  obtained  using  the  absorb¬ 
ing  boundary  condition  to  those  obtained  using  a  standard^ 
NEMD  shock  simulation.  The  only  difference  between  the 
simulations  is  that  the  sample  length  in  the  latter  case  is 
twice  as  long  as  in  the  former.  Both  sets  of  results  corre¬ 
spond  to  t=  1 0.0  ps,  which  is  just  prior  to  the  point  of  maxi¬ 
mum  compression  in  the  longer  cell  (black  lines)  and  4.43  ps 
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EIG.  5.  (Color)  Spatial  profiles  of  (a)  density,  (b)  local  velocity, 
(c)  internal  energy,  and  (d)  temperature  along  the  shock  direction  at 
time  /=  10.0  ps  for  two  simulations  of  shocked  methane  crystal  dif¬ 
fering  only  in  the  initial  sample  length.  Red  lines  correspond  to  the 
simulation  discussed  in  connection  with  Eig.  4,  for  which  the  ab¬ 
sorbing  boundary  condition  was  applied  at  t=5.51  ps;  black  lines 
correspond  to  a  system  twice  as  long  in  the  shock  direction,  which 
has  not  reached  maximum  compression  by  the  end  of  the  simula¬ 
tion.  The  piston  velocity  is  3  km/s  and  the  initial  temperature  is 
50  K. 


after  piston-2  was  applied  in  the  smaller  one  (red  lines). 
Thus,  this  comparison  provides  a  direct  test  of  the  absorbing 
boundary  condition  approach.  While  there  are  observable  de¬ 
viations  for  density,  local  mass  velocity,  internal  energy,  and 
temperature  in  the  regions  of  Fig.  5  for  which  computational 
domains  overlap,  they  are  small  and  insignificant  considering 
their  magnitude  in  light  of  the  pre-  and  postshock  states  of 
the  material. 

One  can  estimate  the  length  of  system  that  would  be  re¬ 
quired,  using  a  standard  NEMD  simulation,  to  obtain  results 
comparable  to  those  shown  in  Eig.  4;  that  is,  one  in  which  a 
shocked  state  is  sustained  for  30  ps.  Eor  Up=3.0  km/s, 
—5.5  ps  is  required  for  the  shock  wave  to  traverse  the 
—47.5  nm  sample  length,  which  means  that  a  simulation  cell 
of  length  —260  nm  would  be  needed.  Simulation  of  a  system 
of  this  size  using  a  complicated  reactive  potential-energy 
function  is  impractical  with  current  computing  capabilities 
even  for  a  single  shock-passage  time,  let  alone  for  the  long 
times  required  to  approach  chemical  equilibrium  (e.g., 
— 150  ps  in  the  following  example,  and  in  many  cases,  much 
longer). 


C.  Chemically  reactive  waves  in  methane 

Analyses  of  shock-induced  chemical  transformations 
were  performed  for  the  cases  Up=^.3  and  1 1. 0  km/s  for 
shocked  liquid  and  solid  methane,  respectively,  based  on  an 
ad  hoc  geometric  definition  of  molecular  connectivity.  Spe¬ 
cifically,  it  was  assumed  that  carbon  atoms  are  chemically 
bonded  when  their  separation  is  within  the  cutoff  distance  for 
the  intramolecular  interactions  in  the  AIREBO  potential 
function.  Thus,  a  molecule  is  defined  as  the  set  of  carbon 
atoms  for  which  any  two  members  of  that  set  can  be  linked 
to  all  other  members  through  an  unbroken  sequence  of 
bonds.  Within  this  framework,  we  define  isolated  carbon 
monomers,  dimers,  trimers,  etc.,  as  having  molecular  sizes  I, 
2,  3,  etc. 

Examination  of  molecular  sizes  after  shock  wave  propa¬ 
gation  through  liquid  methane  for  Up=^.3  km/s  indicates 
that  initial  ethane  production  occurs  with  a  latency  of  —3  ps 
after  shock  wave  passage.  (Recall  that  the  shock  wave  tra¬ 
versal  time  for  the  entire  sample  is  only  about  5.5  ps.)  The 
study  of  several  specific  reaction  events  reveals  a  propensity 
for  ethane  formation  to  occur  by  two  unimolecular  dissocia¬ 
tion  events:  2CH4^2CH3-i-2H  followed  by  recombination 
to  yield  C2H6-fH2.  This  result  is  in  agreement  with  tight- 
binding  MD  simulations^^  and  experimental  measurements 
of  temperature^^  and  conductivity^^  behind  the  shock  front. 
Given  that  molecular  dissociation  is  thermally  activated,  it  is 
not  surprising  that  moderate  differences  between  the  calcu¬ 
lated  and  experimental  pressures  have  comparatively  small 
effects  on  the  predicted  reaction  thresholds  and  chemistry. 

The  decomposition  threshold  for  solid  methane  using 
AIREBO  is  Up=9  km/s,  for  which  the  shock  temperature  is 
—4600  K.  This  temperature  is  close  to  the  value  observed  in 
liquid  methane  shock  propagation,  4400  K  [see  Eig.  4(a)]. 
Elert  et  al?^  carried  out  MD  simulations  of  shock  wave 
propagation  in  solid  methane  and  reported  a  somewhat 
higher  value,  10  km/s,  for  the  decomposition  threshold  ve- 
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FIG.  6.  (Color  online)  Time  dependence  of  molecular  carbon 
cluster  size  in  solid  methane  shocked  with  a  piston  velocity  Up 
=  11  km/s.  The  initial  temperature  and  density  were  50  K  and 
0.53  g/cm^,  respectively.  Curves  C2,  C3,  C4,  and  C5  correspond 
nominally  to  ethane,  propane,  isomers  of  butane,  and  isomers  of 
pentane,  respectively.  The  vertical  dashed  line  at  t=22  ps  indicates 
the  time  at  which  the  second  piston  began  to  move. 

locity.  This  discrepancy  is  likely  caused  by  the  different  ap¬ 
proach  used  for  shock  wave  initiation  in  those  simulations:  a 
finite  flyer  plate  of  several  unit  cells  thickness  was  used, 
which  resulted  in  rarefaction  waves  entering  the  compressed 
region  and,  thus,  decreasing  the  time  available  for  reaction  to 
occur  in  the  fully  compressed  state,  whereas  the  initiating 
piston  used  here  emulates  a  macroscopic  striker. 

In  Fig.  6,  we  show  the  time  evolution  of  carbon  molecular 
sizes  for  solid  methane  shocked  at  Up=\  \  km/s.  The  shock 
wave  traversal  time  is  only  about  2.2  ps  (denoted  by  the 
vertical  line  in  the  figure).  Though  most  of  the  carbon  atoms 
in  the  system  are  members  of  three-atom  molecules  (i.e., 
propanelike  chains,  C3),  about  20%  are  involved  in  clusters 
containing  four  or  more  atoms  and  about  25%  of  all  carbon 
atoms  belong  to  clusters  consisting  of  more  than  50  carbon 
atoms.  Diamond  anvil  cell  experiments  in  which  infrared  ab¬ 
sorption,  Raman  spectroscopy,  and  x-ray  diffraction  were 
measured  indicate  that  diamond  and  hydrogen,  as  well  as 
hydrocarbon  polymer  chains,  are  formed  from  methane  upon 
static  compression  in  the  interval  of  10-50  GPa  and  laser 
heating  to  2000-3000  It  is  likely  that  the  system  stud¬ 
ied  here  at  r=200  ps  represents  an  early  stage  on  the  trans¬ 
formation  path  toward  diamond +H2,  an  overall  process  char¬ 
acterized  by  diffusion-limited  rates  with  time  constants  and 
spatial  scales  that  exceed  those  of  the  present  simulation. 

IV.  SUMMARY  AND  CONCLUSIONS 

We  have  developed  a  practical  approach  for  nonequilib¬ 
rium  molecular  dynamics  simulations  of  shock  waves  that 
allows  the  study  of  shocked  states  for  time  scales  far  larger 
than  the  shock  wave  traversal  time  of  the  MD  simulation 


cell.  A  shock  wave  is  generated  by  driving  a  rigid  piston  into 
a  sample  at  constant  piston  velocity.  A  second  rigid  piston  is 
located  at  the  opposite  end  of  the  simulation  cell,  contiguous 
to  and  equilibrated  with  the  material.  When  the  shock  wave 
reaches  the  end  of  the  simulation  cell,  the  second  piston  be¬ 
gins  to  move  at  the  same  velocity  as  the  first  and,  thus, 
provides  an  “impedance  match”  for  the  shock  wave  across 
the  sample-piston  interface.  With  both  pistons  moving  at  the 
same  velocity,  the  Hugoniot  state  of  the  shocked  material  is 
sustained,  while  the  sample  continues  to  evolve.  This  allows 
a  significantly  longer  simulation  without  significant  interfer¬ 
ence  from  reflected  waves  that  arise  at  interfaces.  The  prin¬ 
cipal  distinction  between  the  present  method  and  the  one  due 
to  Zhao  et  al  is  the  extent  to  which  the  present  one  should 
be  extendable  to  treat  split  wave  structures,  for  instance,  an 
elastic  precursor  followed  by  a  plastic  wave,  or  to  accommo¬ 
date  inhomogeneous  wave  profiles  that  might  require  finite 
piston  acceleration  profiles,  which  themselves  might  vary 
across  the  surface  of  the  piston  (e.g.,  a  shock  wave  propagat¬ 
ing  along  the  longitudinal  axis  of  a  hexagonal-cylindrical 
microphase  segregated  copolymer  morphology). 

We  applied  this  absorbing  boundary  condition  approach  to 
shock  waves  in  methane,  modeled  using  the  AIREBO  force 
field.  We  demonstrated  that  the  method  does  not  introduce 
significant  fluctuations  in  density  or  local  velocity  across  the 
simulation  cell  or  at  the  piston-sample  interfaces,  and  that 
the  internal  energy  (which  is  the  criterion  upon  which  the 
time  to  start  the  second  piston  moving  is  based)  is  equili¬ 
brated  on  a  several  picosecond  time  scale.  We  illustrated  the 
practical  usefulness  of  the  method  by  simulating  shock- 
induced  chemistry  in  methane  on  a  200  ps  time  scale,  for  a 
simulation  cell  with  a  shock  wave  traversal  time  of  only 
2.2  ps.  While  the  simple  implementation  of  the  absorbing 
boundary  condition  described  here  involves  infinite  accelera¬ 
tion  of  the  second  piston,  adaptations  to  provide  more  so¬ 
phisticated  “soft  catches”  should  be  straightforward.  We  ex¬ 
pect  this  general  approach  to  be  useful  for  simulations  of 
shock-induced  dynamics  including  chemistry  in  various  ma¬ 
terials,  with  numerical  integration  stability  and,  of  course, 
the  validity  of  the  classical  approximation  being  the  limiting 
factors. 
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